Exercise

Plot a histogram of the death counts.

hist(GoT$Count)

Because we have count data, it is natural to model the data with a Poisson distribution. If we take the empirical mean of the data as our estimate for \(\lambda\), what does the sampling distribution look like?

y <- GoT$Count
n <- length(y)

# Finish: obtain empirical mean
ybar <- mean(y)
sim_dat <- rpois(n, ybar)
qplot(sim_dat, bins = 20, xlab = "Simulated number of deaths", fill = I("#9ecae1"))

# create data frame for side-by-side histograms
df <- rbind(data.frame(y, "observed") %>% rename(count = 1, data = 2),
            data.frame(sim_dat, "simulated") %>% rename(count = 1, data = 2))
ggplot(df, aes(x = count, fill = data))+
  geom_histogram(position = "dodge", bins = 20) +
  scale_fill_brewer()  +
  labs(x = "Number of deaths", y = "Count")

Plotting the histogram of the data next to the histogram of the sampled data using \(\hat{\lambda} = \bar{y}\), we see that there are many more 0-valued observations in the observed data than there are in the simulated data. We might suspect that the Poisson model will not be a good fit for the data. For now, let’s proceed using the Poisson model and see how it behaves.

Exercise

Use the function mcmc_areas() to plot the smoothed posterior distribution for \(\lambda\) with a 90% Highest Posterior Density region. Changing the prob parameter in mcmc_areas() will change the amount of probability mass that is highlighted. The highlighted region will begin in regions of highest density, and will move towards regions of lower density as prob gets larger.

mcmc_areas(lambda_draws, prob = 0.9)

The line in the plot represents the posterior mean, and the shaded area represents 90% of the posterior density. How does the posterior mean for \(\lambda\) compare to the sample mean?

Exercise

Generate posterior predictive samples using the posterior values of \(\lambda\) and store the result in a length(lambda_draws) by n matrix.

y_rep <- apply(lambda_draws, 1, function(x){rpois(n = n, lambda = x)}) %>% t()

We now step through some graphical posterior predictive checks. Here we plot the histogram of the observed counts against several of the posterior predictions. What do you notice?

ppc_hist(y, y_rep[1:8, ], binwidth = 1)

We can compare the density estimate of \(y\) to the density estimates for several (60) of the \(y_rep\)s. What do you notice here?

ppc_dens_overlay(y, y_rep[1:60, ])

In our first histogram of the data, we noticed the high number of 0 counts. Let us calculate the proportion of zeros in \(y\), and compare that to the propotion of zeros in all of the posterior predictive samples.

prop_zero <- function(x){
  mean(x == 0)
} 
prop_zero(y)

## [1] 0.1643836

ppc_stat(y, y_rep, stat = "prop_zero")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

We can plot the means and variances for all of the simulated datasets, along with the observed statistics.

ppc_stat_2d(y, y_rep, stat = c("mean", "var"))

ppc_error_hist(y, y_rep[1:4, ], binwidth = 1) + xlim(-15, 15)

## Warning: Removed 8 rows containing missing values (geom_bar).

Exercise

Based on on these PPCs, does this model appear to be a good fit for the data?

Exercise

Using the code provided for the simple Poisson model, simluate draws from the posterior density of \(\lambda\) with the “Poisson Hurdle” model. The Stan file you’ll need to use is called poisson-hurdle.stan. Store the resulting object in a variable called fit2.

stan_dat <- list(y = y, N = n)
fit2 <- stan("poisson-hurdle.stan", data = stan_dat, refresh = 0, chains = 2)
lambda_draws2 <- as.matrix(fit, pars = "lambda")

Here we compare the posterior densities of \(\lambda\) from the simple Poisson model and the Poisson hurdle model.

# Extract the sampled vlaues for lambda, and store them in a variable called lambda_draws2:
lambda_draws2 <- as.matrix(fit2, pars = "lambda")

# Compare
lambdas <- cbind(lambda_fit1 = lambda_draws[, 1],
                 lambda_fit2 = lambda_draws2[, 1])

# Shade 90% interval
mcmc_areas(lambdas, prob = 0.9) 

Obtaining posterior samples from the hurdle model is more complicated, so the Stan file includes code to draw from the posterior predictive distribution for you. We extract them here, and store the predictive samples in `y_rep2’.

y_rep2 <- as.matrix(fit2, pars = "y_rep")

Exercise

Use the code given above to produce the same PPC visualizations as before for the new results. Comment on how this second model compares to both the observed data and to the simple Poisson model.

ppc_hist(y, y_rep2[1:8,], binwidth = 1)

ppc_dens_overlay(y, y_rep2[1:60, ])

prop_zero <- function(x){
  mean(x == 0)
} 
ppc_stat(y, y_rep2, stat = "prop_zero")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ppc_stat_2d(y, y_rep2, stat = c("mean", "var"))

ppc_error_hist(y, y_rep2[1:4, ], binwidth = 1) + xlim(-15, 15)

## Warning: Removed 8 rows containing missing values (geom_bar).